
#include "stdafx.hpp"
#include "bitmap.hpp"
#include "grid.hpp"
#include "settings.hpp"
#include "simulator.hpp"
#include "stopwatch.hpp"
#ifdef CAVITY_MULTITHREADING_AVAILABLE
    #include <boost/thread.hpp>
    #include <boost/thread/barrier.hpp>
#endif

void primary_logic();
//void do_relax(vec_of_grids & modes, grid::region & reg, grid & lap);
void analyze(grid & g, grid::region & reg, grid & lap);
void do_freq(grid & g, grid::region & reg);
void do_check(grid & g);
void do_ffactor(grid & g, grid::region & reg);
void do_qfactor(grid & g, grid::region & reg, grid & lap);

int main(int argc, char* argv[])
{
    using namespace settings;
    assert(sizeof(uword) == 4);
    try
    {
		output << "\nHello, physicist. It is now " << DateAndTime() << '\n'
			   << "This is the ADMX-cavity bitmap relaxer. Built at " << (__TIMESTAMP__) << std::endl;
		output << "Command line: ";
		for(unsigned i = 0; i < argc; ++i)
			output << argv[i] << ' ';
		output << '\n' << std::endl;
        set_options(argc, argv);
        check_parameters();
        primary_logic();
		output << "Goodbye, physicist." << std::endl;
    }
    catch(const std::string & s)
    {
        output << "Exception thrown\n"
		  << "\tError: " << s
		  << "\nTerminating..." << std::endl;
    }
    cleanup();

	if(settings::wait)
	{
		output << "Press enter to continue." << std::endl;
		std::cin.get();
	}
}

void primary_logic()
{
    using namespace settings;
    bitmap* bmp = new bitmap(bitmap_path);
    if( !(*bmp) ) throw std::string("Couldn't open bitmap");
    grid::region reg = bmp->region(paths[BOUNDARY]);
    grid::boundary bnd = TEmode ? bmp->boundary() : grid::boundary();
    vec_of_grids modes;
    const uword bmp_height = bmp->height();
    const uword bmp_width = bmp->width();
    delete bmp; //to save memory on small systems with big bitmaps.

    thread_method* tmeth;
    mode_method* mmeth;
 #ifdef CAVITY_MULTITHREADING_AVAILABLE
    if(nthreads>1)
        tmeth = new multi_thread(nthreads,reg,bnd);
    else
 #endif
        tmeth = new single_thread;
    if(mode == 0 && !TEmode)
        mmeth = new low_TM(*tmeth);
    else if(mode == 0 && TEmode)
        mmeth = new low_TE(*tmeth);
    else if(mode > 0 && !TEmode)
        mmeth = new high_TM(*tmeth, modes);
    else if(mode > 0 && TEmode)
        mmeth = new high_TE(*tmeth, modes);
    else
        throw std::string("logic error in primary_logic() -- could not determine mode method.");

    grid lap(bmp_width, bmp_height);     //laplacian field (a back-buffer)
    stopwatch timer(stopwatch::auto_start);
    for(unsigned m = 0; m <= mode; ++m)
    {
	    output << "Beginning mode " << m << ": " << DateAndTime();
 #ifdef CAVITY_MODES_VECTOR_USE_PTR_VECTOR
        modes.push_back( new grid(bmp_width,bmp_height) );
 #else
        modes.push_back( grid(bmp_width,bmp_height) );
 #endif
        grid & g = modes.back();
        //if we're still loading previously calculated modes, do that
        if(static_cast<int>(m) <= settings::load_mode)
        {
              load(g,paths[RESULT],grid::text);
              mmeth->update_sumofsq(g,reg);
              continue;
        }
		settings::change_mode_path(m);  //i.e. "result{n}.dat" -> "result{n+1}.dat"
        functions::set_to setf(0.0);
        special_for_each(g,setf);
        functions::make_rand init_cond(g.height()*g.width());
        g.iterate_region(reg, init_cond);

	 	//if one thread then just run similation here (in parent thread).
	 	if(mode==0 && !TEmode && nthreads <= 1)
	 	{
	 	    std::cout << "Using skipBS mode." << std::endl;
            simulator(*tmeth,*mmeth,g,lap,reg,bnd,settings::n,settings::dt).low_TM_single_thread();
	 	}
        else if(nthreads == 1 || nthreads == 0)
            simulator(*tmeth,*mmeth,g,lap,reg,bnd,settings::n,settings::dt)();
        else //otherwise, do multithread...
        {
          #ifdef CAVITY_MULTITHREADING_AVAILABLE
            boost::thread_group tg;
            for(unsigned i = 0; i < nthreads; ++i)
                tg.create_thread( simulator(*tmeth,*mmeth,g,lap,tmeth->region(i),tmeth->boundary(i),settings::n,settings::dt) );
			tg.join_all();
          #else
            throw std::string("This version was not built for multithreading.");
          #endif
        }

        write_file(g, paths[RESULT].c_str(), grid::text, reg);

        mmeth->update_sumofsq(g,reg);
        analyze(g, reg, lap);
    }
    delete tmeth;
    delete mmeth;
    timer.stop();
    output << DateAndTime() << "That took " << timer.read() << " ticks." << std::endl;
}

void analyze(grid & g, grid::region & reg, grid & lap)
{
    using namespace settings;
    if(check) do_check(g);
    if(frequency || multirun_frequency) do_freq(g,reg);
    if(ffactor) do_ffactor(g,reg);
    if(qfactor) do_qfactor(g,reg, lap);
}

void do_freq(grid & g, grid::region & reg)
{
	//TODO: have do_freq() know about TE modes so that
	// it can maintain the TE boundary conditions. Right
	// now it just relaxes everything, so after a few
	// cycles the TE frequencies are off. The deviation after
	// the first relaxation, however, is negligible.
	using namespace settings;
	std::ofstream out(paths[FREQUENCY].c_str());
	functions::freq freak(out,num_points(reg));
	g.iterate_region(reg,freak);	//calculates frequency at each point
	out.close();

	if(multirun_frequency)
	{
		std::fstream f(paths[MR_FREQ].c_str(),std::ios::out|std::ios::in|std::ios::app);
		if(!f)
            output << "Failed initializing fstream in do_freq()." << std::endl;
		f << "## dimensions of grid used in frequency below: "
		  << std::dec << g.width() << 'x' << g.height() << std::endl;
        f << std::scientific << std::setprecision(grid::precision);

		if(mrindex.empty())	//output: [line number] [frequency]
		{
			f.seekg(0, std::ios::beg);
			unsigned int line = 0;
			while(f.ignore(std::numeric_limits<std::streamsize>::max(), '\n')) ++line;
			f.clear();
			--line;	//we want to start at 0.
			f << line << ' ' << freak.average() << std::endl;
		}
		else	//output: [mrindex] [frequency]
		{
			std::string::size_type index;
			while( (index = mrindex.find(';')) != std::string::npos ) mrindex[index] = ' '; //make every ';' a space
			f << mrindex << ' ' << freak.average() << std::endl;
		}

		f.close();
	}
	output << "Average frequency: " << freak.average() << std::endl;
}

void do_check(grid & g)
{
	//This is some very old code. Checks the centerline of a
	//presumably rectangular configuration against a sine wave.
	using namespace settings;
	std::ofstream out;
	grid::region midline;
	midline.push_back(grid::span(g.width()/2,0,g.height())); //This takes a 1 pixel margin into account.
	//May need to be updated for no margin.
	out.open(paths[CHECK_RATIO].c_str());
	if(extended_check)
	{
		std::ofstream out2(paths[CHECK_ACTUAL].c_str());
		std::ofstream out3(paths[CHECK_EXPECTED].c_str());
		functions::check checkf(out,out2,out3);
		g.iterate_region(midline,checkf);
	}
	else
	{
	    functions::check checkf(out);
	    g.iterate_region(midline,checkf);
	}
	out.close();
}

void do_ffactor(grid & g, grid::region & reg)
{
    functions::sum< identity<grid::data_type> > sumplain;
    g.iterate_region(reg, sumplain);
    printdb( sumplain.value() );

    functions::sum< square<grid::data_type> > sumofsq;
    g.iterate_region(reg, sumofsq);
    printdb( sumofsq.value() );
    output << "Form factor = Square(Sum) / ( Area * Sum(Square) ) = "
	      << ( std::pow(sumplain.value(),grid::data_type(2)) / (num_points(reg) * sumofsq.value()) ) << std::endl;
}

void do_qfactor(grid & g, grid::region & reg, grid & lap)
{
    functions::sum< square<grid::data_type>, true > energy2x_first;
    //Above, the 'true' indicates that the sum be weighted by dielectric constant.
    g.iterate_region(reg, energy2x_first);
    printdb( 2.0*energy2x_first.value() );

    functions::calc_lap lapf(lap);       //function objects
    functions::relax relaxf(lap, settings::dt);
    g.iterate_region(reg, lapf);
    g.iterate_region(reg, relaxf);
    functions::sum< square<grid::data_type>, true > energy2x_second;
    g.iterate_region(reg, energy2x_second);
    printdb( 2.0*energy2x_second.value() );

    //The factors of 2 will cancel in the actual calculation, so I'll omit them
    printdb( 2.0*Pi*energy2x_first.value() / (energy2x_first.value() - energy2x_second.value()) );
    printdb( energy2x_first.value() - energy2x_second.value() );

    std::ofstream out(settings::paths[settings::QFACTOR].c_str());
    if(!out) return;
    unsigned int howmany = 5;
    std::cout << "Calculating quality factor ... ";
    for(unsigned int i = 0; i < howmany; ++i)
    {
        energy2x_first.reset();
        g.iterate_region(reg, energy2x_first);  //get energy
        g.iterate_region(reg, lapf);            //calc laplacian
        g.iterate_region(reg, relaxf);          //apply laplacian
        energy2x_second.reset();
        g.iterate_region(reg, energy2x_second); //get energy again
        out << i << ' ' << ( 2.0*Pi*energy2x_first.value() / (energy2x_first.value() - energy2x_second.value()) ) << std::endl;
    }
    output << "done." << std::endl;
}
